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Summary. Particle suspensions are ubiquitous in our daily life, but are not well 
understood due to their complexity. During the last twenty years, various simulation 
methods have been developed in order to model these systems. Due to varying 
properties of the solved particles and the solvents, one has to choose the simulation 
method properly in order to use the available compute resources most effectively 
with resolving the system as well as needed. Various techniques for the simulation 
of particle suspensions have been implemented at the Institute for Computational 
Physics allowing us to study the properties of clay-like systems, where Brownian 
motion is important, more macroscopic particles like glass spheres or fibers solved 
in liquids, or even the pneumatic transport of powders in pipes. In this paper we will 
present the various methods we applied and developed and discuss their individual 
advantages. 

Key words: Particle suspensions, molecular dynamics, stochastic rotation 
dynamics, lattice Boltzmann method 

1 Introduction 

Adding a fluid to a dry granulate causes the behavior of the mixture to change 
dramatically and a host of unexpected phenomena arises. A good example can 
be studied by anyone on the beach: whereas it is impossible to build a sand 
castle from dry sand, once just a little bit of water has been stirred into the 
sand, one can shape the resulting material almost arbitrarily into surprisingly 
complex arrangements. Adding even more fluid might result in the material 
loosing this stability. If we stir such a mixture, it behaves like a liquid of in- 
creased density Other very common particle-fluid mixtures are ubiquitous in 
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our daily life and include the cacao drink which keeps separating into its con- 
stituents, tooth paste and wall paint which are mixtures of finely ground solid 
ingredients in fluids or blood which is made up of red and white blood cells 
suspended in a solvent. An extreme example is the sand on the beach which 
can be blown away by the wind, ft is important for industrial applications to 
obtain a detailed knowledge of those systems in order to optimize production 
processes or to prevent accidents. 

Long-range fluid-mediated hydrodynamic interactions often dictate the be- 
havior of particle-fluid mixtures. The majority of analytical results for the 
particle scale behavior of suspensions has been obtained in the regime of van- 
ishing Reynolds numbers (viscous flow). For large systems, scientists aim at 
an average, continuum description of the large-scale behavior. However, this 
requires time-consuming and sometimes very difficult experimental measure- 
ments of phenomcnological quantities such as the mean settling speed of a 
suspension, the stress contributions in the system of the individual compo- 
nents (solid and fluid) as functions of, e.g., the solid volume fraction of the 
constituents. 

Computer simulation methods are indispensable for many-particle sys- 
tems, for the inclusion of inertia effects (Reynolds numbers > 1) and Brownian 
motion (Peclet number of order 1). These systems often contain a large num- 
ber of important time scales which differ by many orders of magnitude, but 
nevertheless have to be resolved by the simulation, leading to a large numeri- 
cal effort. However, simulations have the potential to increase our knowledge 
of elementary processes and to enable us to find the aforementioned relations 
from simulations instead of experiments. 

Various simulation methods have been developed to simulate particle-fluid 
mixtures. All of them have their inherent strengths but also some disadvan- 
tages. For example, simplified Brownian Dynamics does not contain long- 
ranged hydrodynamic interactions among particles at all [28]. Brownian Dy- 
namics with full hydrodynamic interactions utilizes a mobility matrix which is 
based on tensor approximations which are exact in the limit of zero Reynolds 
number and zero particle volume fraction [57, 3]. However, the computational 
effort scales with the cube of the particle number due to the inversion of 
matrices. Pair-Drag simulations have been proposed by Silbert et al. [46], 
and include hydrodynamic interactions in an approximative way. They have 
focused on suspensions with high densities (up to 50 %) of uncharged spheri- 
cal colloidal particles. Stokesian Dynamics has been presented by Bossis and 
Brady in the 80s and applied by many authors [21, 1, 69, 34]. For example, 
Melrose and Ball have performed detailed studies of shear thickening colloids 
using Stokesian Dynamics simulations [36, 35]. However, this method is lim- 
ited to Reynolds numbers close to zero and the computational effort is very 
high for dynamical simulations. Even with today's powerful computers it is 
not possible to study the dynamics of more than a few hundred particles. 
The method is still widely used due to its physical motivation and its robust- 
ness, but is complicated to code. Boundary-element methods are more flexible 
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than Stokesian dynamics and can also be used to simulate non-spherical or 
deformable particles, but the computational effort is even higher [48, 43]. 

All these methods assume that hydrodynamic interactions are fully devel- 
oped and that the dynamics of the fluid and of the solved particles can be 
treated as fully separated. In reality, this is not the case. Hydrodynamic inter- 
actions are time dependent due to local stresses at the fluid-particle interfaces. 
A number of more recent methods attempt to describe the time dependent 
long-range hydrodynamics properly with the computational effort scaling lin- 
early with the number of particles. These include recent mesoscopic methods 
like dissipative particle dynamics [14, 13, 7], the lattice-Boltzmann method 
[10, 43, 41, 42, 54, 38], or stochastic rotation dynamics [50, 51, 24, 23]. How- 
ever, for small Reynolds numbers, the computational gain of these methods 
is lost due to the additional effort needed to describe the motion of the fluid. 
Finite element or finite difference methods need a proper meshing of the com- 
putational domain which is not trivial for complicated boundary conditions 
as in the case of dense suspensions. Therefore, many authors only simulated 
a limited number of static configurations rather than the full dynamics of the 
system. Advances in remeshing techniques as well as more powerful computers 
have allowed to overcome these problems. Also, in order to avoid remeshing at 
all, uniform grids can be used [16, 63, 25]. These methods arc flexible and ro- 
bust. They can properly treat non-Newtonian effects and incorporate inertia, 
but are complicated to code. 

For a more detailed description of the simulation methods, experiments or 
theoretical approaches not addressed in this paper, the reader is referred to 
one of the various books on colloid science [49, 44, 64, 53, 62, 26]. 

The remainder of this paper focuses on three different simulation tech- 
niques which have recently been applied to particle-laden flows in our group. 
First, we introduce a method developed by Malevanets and Kapral to model 
a solvent with thermal fluctuations. This approach is used to study the prop- 
erties of claylike colloids [24, 23]. For larger particles, thermal fluctuations are 
undesirable. Here, the lattice Boltzmann method and its extension to particle 
suspension is a very good candidate to study the dynamics of glass spheres 
in a sugar solution [38, 22] . The method is easy to code and has been applied 
to suspensions of spherical and non-spherical particles by various authors. If 
the particles are very massive and the density of the fluid is very low, a full 
hydrodynamic treatment of the solvent is not needed anymore. In the last 
chapter we describe an algorithm based on a coarse-grained description of the 
fluid, so that it is resolved on a length scale larger than the particles. Much 
larger systems can be treated this way, but the coarse-graining is justified only 
in certain situations. As an example we model the pneumatic transport of a 
powder in a pipe which is a common process in many industrial applications 
[52, 66, 67]. 

A more computational demanding and not as easy to code method is a 
Navier Stokes solver for the fluid which is coupled to the particles. The method 
has been successfully applied to the simulation of sedimentation processes 



4 Harting, Hecht, Herrmann, McNamara 

of spherical or non-spherical particles and profits from its well established 
physical background and long standing experience with similar fluid solvers 
in engineering disciplines [72, 63, 25, 40, 17, 18]. These methods have a long 
standing history in our group, but have been described in detail elsewhere and 
will not be covered in this paper. 

2 Simulation of Claylike Colloids: Stochastic Rotation 
Dynamics 

Dense suspensions of small strongly interacting particles are complex systems, 
which are rarely understood on the microscopic level. We investigate proper- 
ties of dense suspensions and sediments of small spherical AI2O3 particles 
by means of a combined Molecular Dynamics (MD) and Stochastic Rotation 
Dynamics (SRD) simulation. Stochastic Rotation Dynamics is a simulation 
method developed by Malevanets and Kapral [50, 51] for a fluctuating fluid. 
The work this chapter is dealing with is presented in more detail in refer- 
ences [24, 23]. 

We simulate clay like colloids, for which in many cases the attractive Van- 
der-Waals forces are relevant. They are often called "peloids" (Greek: clay- 
like). The colloidal particles have diameters in the range of some nm up to 
some Jj.m. The term "peloid"' originally comes from soil mechanics, but parti- 
cles of this size are also important in many engineering processes. Our model 
systems of A^O^-particles of about half a jo.m in diameter suspended in water 
are often used ceramics and play an important role in technical processes. In 
soil mechanics [59] and ceramics science [55] , questions on the shear viscosity 
and compressibility as well as on porosity of the microscopic structure which 
is formed by the particles, arise [73, 47]. In both areas, usually high volume 
fractions (<P > 20%) are of interest. The mechanical properties of these suspen- 
sions are difficult to understand. Apart from the attractive forces, electrostatic 
repulsion strongly determines the properties of the suspension. Depending on 
the surface potential, one can either observe formation of clusters or the parti- 
cles are stabilized in suspension and do sediment only very slowly. The surface 
potential can be adjusted by the pK- value of the solvent. Within Dcbyc-Huckcl 
theory one can derive a so-called 2pK charge regulation model which relates 
the simulation parameters with the pH-value and ionic strength / adjusted 
in the experiment. In addition to the static interactions hydrodynamic effects 
are also important for a complete description of the suspension. Since typical 
Peclet numbers are of order one in our system, Brownian motion cannot be 
neglected. 
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2.1 Molecular Dynamics (MD): Simulation of the Suspended 
Particles 

We study colloidal particles, composing the solid fraction, suspended in a fluid 
solvent. The colloidal particles are simulated with molecular dynamics (MD), 
whereas the solvent is modeled with stochastic rotation dynamics (SRD) as 
described in Sect. 2.2. 

In the MD part of our simulation we include effective electrostatic interac- 
tions and van der Waals attraction, a lubrication force and Hertzian contact 
forces. In order to correctly model the statics and dynamics when approach- 
ing stationary states, realistic potentials are needed. The interaction between 
the particles is described by DLVO theory [28, 61, 47]. If the colloidal par- 
ticles are suspended in a solvent, typically water, ions move into solution, 
whereas their counter ions remain in the particle due to a different resolv- 
ability. Thus, the colloidal particle carries a charge. The ions in solution are 
attracted by the charge on the particles and form the electric double layer. It 
has been shown (see [61]), that the resulting electrostatic interaction between 
two of these particles can be described by an exponentially screened Coulomb 
potential 



where d denotes the particle diameter and r is the distance between the parti- 
cle centers, e is the elementary charge, T the temperature, fee the Boltzmann 
constant, and z is the valency of the ions of added salt. Within DLVO theory 
one assumes linear screening, mainly by one species of ions with valency z 
(e.g. z = +1 for NH4"). The first fraction in (1) is a correction to the original 
DLVO potential, which takes the surface curvature into account and is valid 
for spherical particles [6]. 

The effective surface potential Q is the electrostatic potential at the border 
between the diffuse layer and the compact layer, it may therefore be identified 
with the ^-potential. It includes the effect of the bare charge of the colloidal 
particle itself, as well as the charge of the ions in the Stern layer, where the 
ions are bound permanently to the colloidal particle. In other words, DLVO 
theory uses a renormalized surface charge. This charge can be related to the 
pH value of the solvent within Debye-Hiickel theory [23] . 

£ is the permittivity of the vacuum, e r the relative dielectric constant 
of the solvent, k is the inverse Debye length defined by n 2 = 8tt£bI, with 
the ionic strength / and the Bjerrum length Is- We use 81 for water, which 
implies Ib = 7 A. 

The Coulomb term of the DLVO potential competes with the attractive 
van der Waals term 




(1) 
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An = 4.76 • 10~ 20 J is the Hamaker constant [27] which involves the polar- 
izability of the particles. It is kept constant in our simulations since it only 
depends on the material of the particles and on the solvent. The DLVO po- 
tentials are plotted in Fig. 1 for six typical examples with different depth 
of the secondary minimum. The primary minimum has to be modeled sepa- 
rately, as discussed below. Long range hydrodynamic interactions are taken 




2.6 



(center - center distance) / (particle radius) 



Fig. 1. DLVO potentials for AI2O3 spheres of R = 0.5 um diameter suspended in 
water. These are typical potentials used for our simulations as described below. The 
primary minimum at d/R — 2.0 is not reproduced correctly by the DLVO theory. 
It has to be modeled separately. In most of our cases the existence of the secondary 
minimum determines the properties of the simulated system 



into account in the simulation for the fluid as described below. This can only 
reproduce interactions correctly down to a certain length scale. On shorter 
distances, a lubrication force has to be introduced explicitly in the molecu- 
lar dynamics simulation. The most dominant mode, the so-called squeezing 
mode, is an additional force 

F lub = -(v rel ,f)f 67r7?r ™ d , (3) 

r - r\ - r 2 

with r rod = ; (4) 

ri + r 2 

between two spheres with radii n, r 2 and the relative velocity v re i. rj is the 
dynamic viscosity of the fluid. In contrast to the DLVO potentials the lubri- 
cation force is a dissipative force. When two particles approach each other 
very closely, this force becomes very large. To ensure numerical stability of 
the simulation, one has to limit Fi u b- We do this by introducing a small cutoff 
radius r sc . Instead of calculating Fi u b(f) we take the value for Fi u b(T + r sc ). In 
addition, since the force decays for large particle distances, we can introduce 
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a large cutoff radius r lc for which we assume F lub (r) = if r > r lc . As the 
intention of Fi u b is to correct the finite resolution of the fluid simulation, r sc 
and r\ c have to be adjusted in a way that the dynamic properties, i.e., the 
viscosity of a simulated particle suspension with weak DLVO interactions fit 
the measurements. It turns out that r sc = 1.05(ri + r 2 ) and r\ c = 2.5(ri + r 2 ) 
work best. 

To avoid that the particles penetrate each other, one needs a repulsive force 
depending on their overlap. We are using a Hertz force described by the po- 
tential 



where K could be expressed by the elastic modulus of A1 2 3 . This would 
determine the simulation time step, but to keep the computational effort 
relatively small, we determine the time step using the DLVO-potentials as 
described later on and then choose a value for K. Two aspects have to be 
considered: K has to be big enough so that the particles do not penetrate 
each other by more than approximately 10% and it may not be too big, so 
that numerical errors are kept small, which is the case when the collision time 
is resolved with about 20 time steps. Otherwise total energy and momentum 
are not conserved very well in the collision. 

The Hertz force also contains a damping term in normal direction, 



with a damping constant /3 and for the transverse direction a viscous friction 
proportional to the relative velocity of the particle surfaces is applied. 
Since DLVO theory contains the assumption of linear polarizability, it holds 
only for large distances, i.e., the singularity when the two spheres touch does 
not exist in reality. Nevertheless, there is an energy minimum about 30 k^T 
deep, so that particles which come that close would very rarely become free 
again. To obtain numerical stability of our simulation, we model this mini- 
mum by a parabolic potential, some k^T deep (e.g. SksT). The depth of the 
minimum in our model is much less than in reality, but the probability for 
particles to be trapped in the minimum has to be kept low enough so that 
only few of them might escape during simulation time. 

For the integration of the translational motion we utilize a velocity Verlet 
algorithm [4] to update the velocity and position of particle i according to the 
equations 



Vkcrtz = K(d - rf' 2 if r<d 



(5) 



F D amp = -(v rc i,f)f/Vf ~n-r 2 



(6) 



p 

Xi(i + St) = Xi(t) + 5i Vi(t) + St 2 — 




(7) 



(8) 



For the rotation, a simple Euler algorithm is applied: 
Wi(t + 6t) = uii{t) + 8tTi , 



(9) 
(10) 
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where Wj(t) is the angular velocity of particle i at time t, Tj is the torque 
exerted by non central forces on the particle i, is the orientation of 

particle i at time t, expressed by a quaternion, and F(di,uJi, St) gives the 
evolution of of particle i rotating with the angular velocity u>i (t) at time t. 
The concept of quaternions [4] is often used to calculate rotational motions 
in simulations, because the Euler angles and rotation matrices can easily be 
derived from quaternions. Using Euler angles to describe the orientation would 
give rise to singularities for the two orientations with i? = ±90°. The numerical 
problems related to this fact and the relatively high computational effort of a 
matrix inversion can be avoided using quaternions. 

2.2 Stochastic Rotation Dynamics (SRD): Simulation of the Fluid 

The Stochastic Rotation Dynamics method (SRD) introduced by Malevanets 
and Kapral [50, 51] is a promising tool for a coarse-grained description of a 
fluctuating solvent, in particular for colloidal and polymer suspensions. The 
method is also known as "Real-coded Lattice Gas" [33] or as "multi-particlc- 
collision dynamics" (MPCD) [60]. It can be seen as a "hydrodynamic heat 
bath" , whose details are not fully resolved but which provides the correct 
hydrodynamic interaction among embedded particles [45]. SRD is especially 
well suited for flow problems with Peclet numbers of order one and Reynolds 
numbers on the particle scale between 0.05 and 20 for ensembles of many 
particles. The method is based on so-called fluid particles with continuous 
positions and velocities. Each time step is composed of two simple steps: One 
streaming step and one interaction step. In the streaming step the positions 
of the fluid particles are updated as in the Euler integration scheme known 
from Molecular Dynamics simulations: 

Y l {t + T) =n(t)+TVi{t)i , (11) 

where Ti(t) denotes the position of the particle i at time t, Vj(i) its velocity 
at time t and r is the time step used for the SRD simulation. After updating 
the positions of all fluid particles they interact collectively in an interaction 
step which is constructed to preserve momentum, energy and particle number. 
The fluid particles are sorted into cubic cells of a regular lattice and only the 
particles within the same cell are involved in the interaction step. First, their 
mean velocity Uj(t') = N ^ t ,) Si=a* ^ v i(*) i s calculated, where Uj(t') denotes 
the mean velocity of cell j containing Nj(t') fluid particles at time t' = t + r. 
Then, the velocities of each fluid particle in cell j are updated as: 

v 4 (t + r) = u 3 (t') + n 3 (t') ■ [ Vi (i) - uj(t')] . (12) 

£lj(t') is a rotation matrix, which is independently chosen randomly for each 
time step and each cell. We use rotations about one of the coordinate axes 
by an angle ±a, with a fixed. This has been suggested by M. Straufi in [70]. 
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The coordinate axis as well as the sign of the rotation are chosen by random, 
resulting in six possible rotation matrices. The mean velocity Uj(t) in the cell 
j can be seen as streaming velocity of the fluid at the position of the cell j at 
the time t, whereas the difference [vj(t) — Uj(t')} entering the interaction step 
can be interpreted as a contribution to the thermal fluctuations. 

In order to remove low temperature anomalies and to achieve exact 
Galilean-invariance, we use a modification of the original algorithm [29]: all 
particles are shifted by the same random vector with components in the inter- 
val [—a/2, a/2] before the collision step. Particles are then shifted back by the 
same amount after the collision. The random vectors of consecutive iterations 
are uncorrelated. Ihle and Kroll have discussed in [30, 31] why this simple pro- 
cedure works and shown that it leads to transport coefficients independent of 
an imposed homogeneous flow field. In [32] and [37] analytical calculations of 
the transport coefficient of this method are presented. 

Two different methods to couple the SRD and the MD simulation have 
been introduced in the literature. Inoue et al. proposed a way to implement no 
slip boundary conditions on the particle surface [33], whereas Falck et al. [15] 
have developed a "more coarse grained" method we describe shortly in the 
following section. 

2.3 Coupling of the MD and the SRD Simulation Part 

To couple the two parts of the simulation, MD on the one hand and SRD on 
the other one, the colloidal particles are sorted into the SRD boxes and their 
velocities are included in the rotation step. This technique has been used to 
model protein chains suspended in a liquid [15, 74]. Since the mass of the fluid 
particles is much smaller than the mass of the colloidal particles, one has to 
use the mass of each particle — colloidal or fluid particle — as a weight factor 
when calculating the mean velocity 

Uj(t') = TTTTTT E v '( ( ) m > ' ( 13 ) 
^ ' »=1 

with Mj(t') = m i > ( 14 ) 

i=i 

where we sum over all colloidal and fluid particles in the cell, so that Nj(t') 
is the total number of both particles together, rrik is the mass of the particle 
with index i and therefore Mj(t') gives the total mass contained in cell j at 
the time t' = t + r. 
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Fig. 2. Schematic phase diagram for volume fraction $ = 35% in terms of pH-value 
and ionic strength involving three different phases: a clustering regime due to van 
der Waals attraction, stable suspensions where the charge of the colloidal particles 
prevents clustering, and a repulsive structure for further increased electrostatic re- 
pulsion. This work concentrates on state A (pH = 6, I — 3mmol/l) in the suspended 
phase, state B (pH = 6, I = 7mmol/l) close to the phase border but already in the 
clustered phase, and state C (pH = 6, I = 25mmol/l) in the clustered phase [23] 

2.4 Results 
Phase Diagram 

Depending on the experimental conditions, one can obtain three different 
phases: A clustered region, a suspended phase, and a repulsive structure. 
These phases can be reproduced in the simulations and we can quantitatively 
relate interaction potentials to certain experimental conditions. A schematic 
picture of the phase diagram is shown in Fig. 2. Close to the isoelectric point 
(pR = 8.7), the particles form clusters for all ionic strengths since they are 
not charged. At lower or higher pR values one can prepare a stable suspension 
for low ionic strengths because of the charge, which is carried by the colloidal 
particles. At even more extreme pR values, one can obtain a repulsive struc- 
ture due to very strong electrostatic potentials (up to £ = 170 mV for pR = 4 
and / = lmmol/1, according to our model). The repulsive structure is char- 
acterized by an increased shear viscosity. In the following we focus on three 
states: State A (pR = 6, I = 3mmol/l) is in the suspended phase, state B 
(pR = 6, 1 = 7mmol/l) is a point already in the clustered phase but still close 
to the phase border, and state C (pR = 6, I = 25mmol/l) is located well in 
the clustered phase. 

Some typical examples for the different phases are shown in Figs. 3a)-d). 
These examples are meant to be only illustrative and do not correspond ex- 
actly to the cases A-C in Fig. 2 denoted by uppercase letters. In the suspended 
case (a), the particles are mainly coupled by hydrodynamic interactions. One 
can find a linear velocity profile and a slight shear thinning. If one increases 
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c) central cluster d) plug flow 

Fig. 3. Images of four different cases. For better visibility we have chosen smaller 
systems than we usually use for the calculation of the viscosity. The colors denote 
velocities: Dark particles are slow, bright ones move fast. The potentials do not 
correspond exactly to the cases A-C in Fig. 2, but they show qualitatively the 
differences between the different states: a) Suspension like in state A, at low shear 
rates, b) Layer formation, which occurs in the repulsive regime, but also in the 
suspension (state A) at high shear rates, c) Strong clustering, like in state C, so 
that the single cluster in the simulation is deformed, d) Weak clustering close to 
the phase border like in state B, where the cluster can be broken into pieces, which 
follow the flow of the fluid (plug flow) 

the shear rate 7 > 500/s, the particles arrange in layers. The same can be ob- 
served if the Debye-screening length of the electrostatic potential is increased 

(b) , which means that the solvent contains less ions (/ < 0.3 mmol/1) to screen 
the particle charges. On the other hand, if one increases the salt concentration, 
electrostatic repulsion is screened even more and attractive van der Waals in- 
teraction becomes dominant (/ > 4 mmol/1). Then the particles form clusters 

(c) , and viscosity rises. A special case, called "plug flow", can be observed for 
high shear rates, where it is possible to tear the clusters apart and smaller 
parts of them follow with the flow of the solvent (d). This happens in our sim- 
ulations for / = 25 mmol/1 (state C) at a shear rate of 7 > 500/s. However, 
as long as there are only one or two big clusters in the system, it is too small 
to expect quantitative agreement with experiments. In these cases we have to 
focus on state B (I = 7 mmol/1) close to the phase border. 

We restrict ourselves to the region around pH = 6 where we find the 
phase border between the suspended region and the clustered regime at about 
I = 4 mmol/1 in the simulations as well as in the experiments. Also the shear 
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rate dependence of the viscosity is comparable in simulations and experiments 
as discussed below. 

Shear Profile and Shear Viscosity 

In each of the three phases a typical velocity profile of the shear flow occurs. In 
the suspended phase one finds a linear velocity profile (Fig. 4a)) with nearly 
Newtonian flow. The particles are distributed homogeneously thus the density 
profile is structureless (Fig. 5a)). The motion of the particles is only weakly 
coupled by the hydrodynamic forces. At high enough shear rates (7 > 500) 
the particles arrange in layers parallel to the shear plane, as can be seen in the 
density profile Fig. 5b), too. This arrangement minimizes collisions between 
the particles. As a result, the shear viscosity descents as shown in Fig. 6, 
which we discuss in more detail below. Shear induced layer formation has 
been reported in literature for different systems. Vermant and Solomon have 
reviewed this topic recently [71]. 
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Fig. 4. Profiles of tangential velocity component in normal direction: a) Linear 
profile in the suspended regime, state A of Fig. 2 (I — 3mmol/l) at 7 = 500/s) 

b) Cluster formation in state C (I = 25mmol/l) at 7 = 100/s. In principle one 
could determine the viscosity of one single cluster from the central plateau, but this 
is not the viscosity found in experiments. There, one measures the viscosity of a 
paste consisting of many of these clusters 

c) Same as case b) but with higher shear rate (500/s). Hydrodynamic forces are 
large enough to break the cluster into two pieces. The velocity axis is scaled with 
the shear velocity vs for better comparability 

In the clustered phase, the clusters move in the fluid as a whole. They 
are deformed, but since the inter-particle forces are stronger than the hydro- 
dynamic forces, the cluster moves more like a solid body than like a fluid. 
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Fig. 5. Density profiles: a) Suspended case: State A in Fig. 2 (I — 3mmol/l), at 
low shear rates (7 = 50/s). The density distribution is homogeneous 

b) Shear induced layer formation: This is state A as in graph a) of this figure, but 
for a high shear rate (7 = 1000/s) 

c) Strong attractive forces in state C (I = 25 mmol/1): For low shear rates (7 = 50/s) 
only one central cluster is formed, which is deformed slowly 
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Fig. 6. Comparison between simulation and experiment: viscosity in dependence of 
the shear rate for the states A (I = 3 mmol/l) and B (I — 7 mmol/l) of Fig. 2. Note: 
shear thinning is more pronounced for the slightly attractive interactions in state 
B than for the suspended state A. Lines denote experimental data [58], points are 
results from our simulations 



Often there is one big cluster that spans the whole system. The density pro- 
file (Fig. 5c)) increases in the central region and decays at the regions close to 
the border, since particles from there join the central cluster. When averaging 
the velocity profile in the shear flow, one finds a very small velocity gradient 
in the center of the shear cell and fast moving particles close to the wall, 
where the shear is imposed (Fig. 4b)). The velocity profile is non-linear on 
the length scale of the simulations. In the experiment the physical dimensions 
are much larger and therefore the velocity profile can become approximately 
linear again if the system consists of many large clusters. However, due to 
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the computational effort in simulations it is today impossible to measure the 
shear viscosity for these strongly inhomogeneous systems. 

Closer to the phase border clusters can then be broken up into small pieces 
by the hydrodynamic forces at least for high shear rates. In state C of Fig. 2 
this happens for the first time at 7 — 500/s, so that one can find two clusters 
in the system moving in opposite directions. The velocity profile of this case 
is shown in Fig. 4c). For even higher shear rates or closer to the phase border 
(e.g. state B), the clusters are broken into smaller pieces. Then, they move 
in the shear flow with an approximately linear velocity profile. Due to van 
der Waals attraction the system resists with stronger shear forces and the 
viscosity is higher than in the suspended case (Fig. 6). 

In Fig. 6 the simulation results are shown together with the experimen- 
tal results, both for the two cases of a slightly clustered system in state B 
(7 = 7mmol/l) and a suspension (state A, I = 3mmol/l). For the suspension 
(state A) the viscosity decreases with the shear rate ("shear thinning"). The 
experimental data and the simulation are consistent within the accuracy of 
our model. There are several reasons for which our model does not fit exactly 
the measurements: Even though we use a charge regulation model to deter- 
mine the input parameters for the DLVO potentials, microscopic properties 
like the surface density of sites, where ions can be adsorbed on the surface of 
the colloidal particle, have to be determined indirectly. Measurements of the 
(■-potential in certain conditions provide data to fit the unknown microscopic 
parameters. Furthermore, we have monodisperse spheres, which is another 
simplification in our model. 

For the slightly clustered case (state B) an increase of the shear viscos- 
ity, compared to the suspended case, can be observed in the experiment as 
well as in the simulations. Shear thinning becomes more pronounced, because 
clusters are broken up, as mentioned above. However, the shear rate depen- 
dence is stronger in the simulations than in the experiment. This can be the 
first indication of finite size effects. We have studied the dependence of the 
simulated shear viscosity in dependence of the system size. The effect is most 
important for low shear rates. 

3 Transport Phenomena and Structuring in Suspensions: 
Lattice-Boltzmann Simulations 

For industrial applications, systems with rigid boundaries, e.g. a pipe wall, are 
of particular interest since structuring effects might occur in the solid fraction 
of the suspension. Such effects are known from dry granular media resting on 
a plane surface or gliding down an inclined chute [56, 68]. In addition, the 
wall causes a demixing of the solid and fluid components which might have 
an unwanted influence on the properties of the suspension. Near the wall one 
finds a thin lubrication layer which contains almost no particles and causes a 
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so-called "pseudo wall slip" . Due to this slip the suspension can be transported 
substantially faster and less energy is dissipated. 

We expect structuring close to a rigid wall at much smaller concentrations 
than in granular media because of long-range hydrodynamic interactions. In 
[38], we study these effects by the means of particle volume concentrations 
versus distance to the wall. 

3.1 The Lattice-Boltzmann Method 

The latticc-Boltzmann method is a simple scheme for simulating the dynamics 
of fluids. By incorporating solid particles into the model fluid and imposing the 
correct boundary condition at the solid/fluid interface, colloidal suspensions 
can be studied. Pioneering work on the development of this method has been 
done by Ladd et al. [41, 42, 43] and we use their approach to model sheared 
suspensions near solid walls. 

The lattice-Boltzmann (hereafter LB) simulation technique which is based 
on the well-established connection between the dynamics of a dilute gas and 
the Navier-Stokes equations [9]. We consider the time evolution of the one- 
particle velocity distribution function n(r, v, t), which defines the density of 
particles with velocity v around the space-time point (r, t). By introducing the 
assumption of molecular chaos, i.e. that successive binary collisions in a dilute 
gas are uncorrelated, Boltzmann was able to derive the integro-diffcrential 
equation for n named after him [9] 



where the left hand side describes the change in n due to collisions. 

The LB technique arose from the realization that only a small set of dis- 
crete velocities is necessary to simulate the Navier-Stokes equations [20] . Much 
of the kinetic theory of dilute gases can be rewritten in a discretized version. 
The time evolution of the distribution functions n is described by a discrete 
analogue of the Boltzmann equation [43] : 



where Ai is a multi-particle collision term. Here, rii(r,t) gives the density 
of particles with velocity Ci at (r,t). In our simulations, we use 19 different 
discrete velocities c,. The hydrodynamic fields, mass density g, momentum 
density j = gu, and momentum flux 77, are moments of this velocity distri- 
bution: 




(15) 



n,(r + aAt, t + At)= m(r, t) + A(r, t) , 



(16) 






(17) 



We use a linear collision operator, 
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AM) =M ij (n j -nf) , 



(18) 



where Mij = az ^^" — - is the collision matrix and the equilibrium distri- 
bution [10], which determines the scattering rate between directions i and j. 
For mass and momentum conservation, My satisfies the constraints 



M 
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E M « = 



E e l M lJ = . 



(19) 



i=l 



We further assume that the local particle distribution relaxes to an equilibrium 
state at a single rate r and obtain the lattice BGK collision term [5] 



(20) 



By employing the Chapman-Enskog expansion [9, 19] it can be shown that 
the equilibrium distribution 
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and the kinematic viscosity [43] 



V 



2t- 1 



Qf 9 

properly recovers the Navier-Stokes equations 
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(21) 
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(24) 



3.2 Fluid-Particle Interactions 



To simulate the hydrodynamic interactions between solid particles in sus- 
pensions, the lattice-Boltzmann model has to be modified to incorporate the 
boundary conditions imposed on the fluid by the solid particles. Stationary 
solid objects are introduced into the model by replacing the usual collision 
rules (Equation (20)) at a specified set of boundary nodes by the "link-bounce- 
back" collision rule [54]. When placed on the lattice, the boundary surface cuts 
some of the links between lattice nodes. The fluid particles moving along these 
links interact with the solid surface at boundary nodes placed halfway along 
the links. Thus, a discrete representation of the surface is obtained, which be- 
comes more and more precise as the surface curvature gets smaller and which 
is exact for surfaces parallel to lattice planes. 
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Since the velocities in the lattice-Boltzmann model are discrete, boundary 
conditions for moving suspended particles cannot be implemented directly. 
Instead, we can modify the density of returning particles in a way that the 
momentum transferred to the solid is the same as in the continuous velocity 
case. This is implemented by introducing an additional term Ab in (16) [41]: 

Aba = ^'T' Cl . (25) 

c s 

with c s being the velocity of sound and coefficients uj Ci from (22). 

To avoid redistributing fluid mass from lattice nodes being covered or un- 
covered by solids, we allow interior fluid within closed surfaces. Its movement 
relaxes to the movement of the solid body on much shorter time scales than 
the characteristic hydrodynamic interaction [41]. 

If two particle surfaces approach each other within one lattice spacing, 
no fluid nodes are available between the solid surfaces. In this case, mass is 
not conserved anymore since boundary updates at each link produce a mass 
transfer Aba 3 (a =cell size) across the solid-fluid interface [41]. The total 
mass transfer for any closed surface is zero, but if some links are cut by two 
surfaces, no solid-fluid interface is available anymore. Instead, the surface of 
each particle is not closed at the solid-solid contacts anymore and mass can 
be transferred in-between suspended particles. Since fluid is constantly added 
or removed from the individual particles, they never reach a steady state. In 
such cases, the usual boundary-node update procedure is not sufficient and a 
symmetrical procedure which takes account of both particles simultaneously 
has to be used [42] . Thus, the boundary-node velocity is taken to be the aver- 
age of that computed from the velocities of each particle. Using this velocity 
the fluid populations are updated (Equation (25)), and the force is computed; 
this force is then divided equally between the two particles. 

If two particles are in near contact, the fluid flow in the gap cannot be 
resolved by LB. For particle sizes used in our simulations (R < 5a), the lubri- 
cation breakdown in the calculation of the hydrodynamic interaction occurs 
at gaps less than 0.1R [54]. This effect "pushes" particles into each other. 

To avoid this force, which should only occur on intermolccular distances, 
we use a lubrication correction method described in [54]. For each pair of 
particles a force 

„ R1R2 ( 1 1\ ria , , .... 

Fiub = -67T7? — — - — I U12 • 1 , , h <h N (26) 



{R l +R 2 ) 2 \h h N J |n 2 

is calculated, where ui 2 = ui — u 2 , h = |ri 2 | — R\ — i? 2 is the gap between the 
two surfaces and a cut off distance Hn — |a [43]. For particle- wall contacts 
we apply the same formula with R 2 — > 00 and h = |ri 2 | — R\. The tangential 
lubrication can also be taken into account, but since it has a weaker logarith- 
mic divergence and its breakdown does not lead to serious problems, we do 
not include it in our simulations. 
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3.3 Particle Motion 

The particle position and velocity are calculated using Newton's equations 
in a similar manner as in section on SRD simulations. To avoid repetition, 
the reader is referred to Sect. 2.1. However, particles do not feel electrostatic 
interactions, but behave like hard spheres in the case presented in this section. 

3.4 Simulations 

The purpose of our simulations is the reproduction of rheological experiments 
on computers. We simulate a representative volume element of the experimen- 
tal setup and compare our calculations with experimentally accessible data, 
i.e. density profiles, time dependence of shear stress and shear rate. We also 
get experimentally inaccessible data from our simulations like translational 
and rotational velocity distributions, particle-particle and particle-wall inter- 
action frequencies. The experimental setup consists of a rheoscope with two 
spherical plates, which distance can be varied. The upper plate can be rotated 
cither by exertion of a constant force or with a constant velocity, while the 
complementary value is measured simultaneously. The material between the 
rheoscope plates consist of glass spheres suspended in a sugar-water solution. 
The radius of the spheres varies between 75 and 150 jim. For our simulations 
we assume an average particle radius of 112.5 |tm. The density and viscosity 
of the sugar solution can also be changed. We simulate only the behavior of 
a representative volume clement which has the experimental separation be- 
tween walls, but a much lower extension in the other two dimensions than the 
experiment. In these directions we employ periodic boundary conditions for 
particles and for the fluid. 

Shearing is implemented using the "link-bounce-back" rule with an addi- 
tional term A^^ at the wall in the same way as already described for particles 
(Equation (25) with now being the velocity of the wall). 

To compare the numerical and experimental results, we need to find char- 
acteristic dimensionless quantities of the experiment which then determine 
the simulation parameters. For this purpose we use the ratio of the rheoscope 
height and the particle size A, the particle Reynolds number 3? and the volume 
fraction of the particles (f>. The simulation results are provided with units by 
calculating the length of the lattice constant a and the duration of one time 
step as described in [38]. 

3.5 Results 

Figure 7 shows a snapshot of a suspension with 50 spheres after 5772500 time 
steps which are equivalent to 729 s. The vector g represents the direction of 
gravity and vj depicts the velocity of the sheared wall. 

The particles feel a gravitational acceleration g = 0.8 m/s 2 , have a mass 
m = 7.7 • 10~ 8 kg, a Reynolds number 3? = 4.066875 • 10~ 4 , and a radius 
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Fig. 7. A snapshot of a suspension with 50 spheres (radius R = 1.125 -10 4 m, 
mass m = 7.7 • 10 -8 kg) at time t = 729 s. The volume of the simulated system is 
1.83 • 1(T 3 x 1.83 • 1(T 3 x 3.375 • 1(T 3 m = 11.3025 • 10" 9 m 3 , acceleration of gravity 
g — 0.80 m/s 2 , and shear velocity v s = 3.375 • 10~ 2 m/s. The fluid has a viscosity 
r\ — 450 mPa • s and density Qf = 1446 v p-. This visualization is a typical example 
for a system that has reached a steady state: All particles have fallen to the ground 
due to the exerted gravitational force and most of the system has no particles [38] 



R = 1.125 • 1(T 4 m. The system size is 1.83 • 10" 3 x 1.83 • 1CT 3 x 3.375 • 1CT 3 m 
which corresponds to a lattice size of 32 x 32 x 59. The density of the fluid is set 
to Qf = 1446^ and its viscosity is rj — 450 mPa • s. The walls at the top and 
the bottom are sheared with a relative velocity v s — 3.375- 10 -2 m/s. Figure 7 
is a representative visualization of our simulation data and demonstrates that 
after the system has reached its steady state, all particles have fallen to the 
ground due to the influence of the gravitational force. Most of the simulation 
volume is free of particles. 

In order to quantitatively characterize structuring effects, we calculate the 
particle density profile of the system by dividing the whole system into layers 
parallel to the walls and calculating a partial volume for each particle 
i crossing such a layer j. The scalar Vij is given by the volume fraction of 
particle i that is part of layer j: 

Vij = tt (i? 2 (R™* R™ n ) - 1 (R™* R^) ) (27) 

If the component perpendicular to the wall of the radius vector of the 
center of sphere i lies between r™ m and r™ ax , we have 

rf n = [j ~^AL Z -R, 
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Fig. 8. Density profiles from simulations with two different shear rates 7 = 10 s _1 
(a) and 7=1 s _1 (b). Other parameters are equal to those given in Fig. 7. (a) shows 
five peaks with separations about one particle diameter, which reveal the forming 
of particle layers. The number of particles per layer is decreasing with increasing 
distance to the wall, and the change in particle numbers is caused by gravity which 
is directed perpendicular to the wall at z — 0. Although we used the same gravity 
and particle numbers, there are only three peaks in (b) and their width is higher 
than in (a), demonstrating that the structuring effects strongly relate to the shear 
rate 



and 
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Finally, the sum of all weights associated with a layer is divided by the volume 
of the layer 



L 2 



" J L x -L y -AL Z ^ M 

with L Xl L y being the system dimensions between periodic boundaries, L z 
the distance between walls, M the number of layers, and AL Z the width of a 
single layer. 

Density profiles calculated by this means for systems with two different 
shear rates 7 = 10 s _1 and 7 = 1 s _1 are presented in Fig. 8. All other 
parameters are equal to the set given in the last paragraph. The peaks in Fig. 
8 demonstrate that at certain distances from the wall the number of particles 
is substantially higher than at other positions. The first peak in both figures is 
slightly below one particle diameter, which can be explained by a lubricating 
fluid film between the first layer and the wall which is slightly thinner than 
one particle radius. Due to the small amount of particles, time dependent 
fluctuations of the width of the lubricating layer cannot be neglected and 
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a calculation of the exact value is not possible. The five peaks in Fig. 8a 
have similar distances which are equal to one particle diameter. These peaks 
can be explained by closely packed parallel layers of particles. Due to the 
linear velocity profile in z-direction of the fluid flow, every layer adopts the 
local velocity of the fluid resulting in a relative velocity difference between 
two layers of about 2RJ. These layers stay stable in time with only a small 
number of particles being able to be exchanged between them. 

Figure 8b only shows three peaks with larger distances than in Fig. 8a. 
However, the average slope of the profile is identical for both shear rates. For 
smaller shear rates, velocity differences between individual layers are smaller, 
too. As a result, particles feel less resistance while moving from one layer to 
another. Every inter-layer transition distorts the well defined peak structure 
of the density distribution resulting in only three clearly visible peaks in Fig. 
8b. 

With changing time, the first peak stays constant for both shear rates. The 
shape, number and position of all other peaks is slightly changing in time. 



t=28000 




Fig. 9. A snapshot of a suspension with 1536 spheres after 28000 timesteps used to 
gain statistics of particle velocity distributions 

We are currently investigating the occurrence of non-Gaussian velocity 
distributions of particles for higher particle densities and higher shear rates. 
For this, improvements of the method are mandatory in order to prevent 
instabilities of the simulation. By utilizing an implicit scheme for the update 
of the particle velocities [43, 54] we are able to overcome artefacts caused by 
numerical inaccuracies at high volume fractions or shear rates. Figure 9 shows 
a snapshot of a system containing 1536 particles after 28000 timesteps. 
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The lattice Boltzmann has been extended in order to include thermal 
fluctuations [43, 2] . With these modifications the method is another candidate 
to simulate suspensions where Brownian motion cannot be neglected. 

4 Plug Conveying in Vertical or Horizontal Tubes: a 
Coarse Grained Model for the Fluid Flow 

4.1 Model Description 

Another approach to modeling two phase flow is to course-grain the fluid, so 
that it is resolved on a length scale larger than the grains. The advantage is 
that much larger systems can be treated, but the disadvantage is that this 
coarse-graining is justified only in certain situations. One of those situations is 
when the density of the fluid is small compared to that of the grains, and the 
Reynolds number of the grains is small. It is then possible to neglect the inertia 
of the fluid, which means that all momentum is contained in the grains. The 
fluid transfers momentum between grains, but stores no momentum itself. And 
when the Reynolds number of the grains is small, one can treat the granulate 
phase as a moving porous medium. In the following, we present the model in 
more detail. 

Gas Model 

The model for the gas simulation was first introduced by McNamara and 
Flekk0y [52] and has been implemented for the two-dimensional case to sim- 
ulate the rising of bubbles within a fluidized bed. We developed a three- 
dimensional version of this algorithm. 

The algorithm is based on the mass conservation of the gas and the granu- 
lar medium. Conservation of grains implies that the density g p of the granular 
medium obeys 

^ + V -(U0 P ) = O, Q P = Q s {l-4>), (29) 

where g p is the mass density of the material making up the particles, the 
porosity of the medium is <f> (i.e. the fraction of the space available to the 
gas), and the velocity of the granulate is u. 

The mass conservation equation for the gas is 

+ V • (v g g g ) = , g g cx , (30) 

where g g is the mass density of the gas averaged over the total volume of the 
granular medium and v g its velocity. This equation can be transformed into 
a differential equation for the gas pressure P using the ideal gas equation, 
together with the assumption of uniform temperature. 
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The velocity v g of the gas is related to the granulate velocity u through 
the d'Arcy relation: 

- VP =^ V *- U >' (31) 

where 77 is the dynamic viscosity of the air and k is the permeability of the 
granular medium. This relation was first given by d'Arcy in 1856 [12]. The 
d'Arcy relation is preferred here over the Ergun equation, because it is linear 
in the velocity. This makes the simplification steps done later possible. For the 
permeability k the Carman-Kozeny relation [8] was chosen, which provides a 
relation between the porosity <f), the particle diameter d and the permeability 
of a granular medium of monodisperse spheres, 

d 2 <j> 3 , s 

= Tso(^W ' (32) 

Combining (29), (30) and (31) results in a nonlinear differential equation for 
the gas pressure: 

^+uVF)=V(P^VF)-PVu. (33) 

After linearizing around the normal atmospheric pressure Po the resulting dif- 
ferential equation only depends on the relative pressure P' (P = Pq + P'), the 
porosity <j) and the granular velocity u, which can be derived from the particle 
simulation, and three constants: the viscosity 77, the particle diameter d and 
the pressure P®. 

r)V" P P 

This differential equation can be interpreted as a diffusion equation with a 
diffusion constant D = (pK((j))/r]. The equation is solved numerically, using a 
Crank-Nickelson approach for the discretization. Each dimension is integrated 
separately. 



Granulate Algorithm 

The model for the granular medium simulates each grain individually us- 
ing a discrete element simulation (DES). For the implementation of the dis- 
crete element simulation we used a version of the molecular dynamics method 
described by Cundall [11]. The particles are approximated as monodisperse 
spheres, rotations in three dimensions are taken into account. 
The equation of motion for an individual particle is 

^ mVP ,„„, 

mx = mg + F c — , (35) 

Qs(l - <j>) 



where m is the mass of a particle, g the gravitation constant and F c the 
sum over all contact forces. The last term, the drag force, is assumed to be a 
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volume force given by the pressure drop VP and the local mass density of the 
granular medium g s (l — cf>), which is valid for monodisperse granular media. 

The interaction between two particles in contact is given by two force 
components: a normal and a tangential component with respect to the particle 
surface. The normal force is the sum of a repulsive clastic force (Hooke's law) 
and a viscous damping. The tangential force opposes the relative tangential 
motion and is proportional to the normal force (sliding Coulomb friction) or 
proportional to the relative tangential velocity (viscous damping). Viscous 
damping is used only for small relative tangential velocities. 

Gas-Grain Interaction 

The simulation method uses both a continuum and a discrete element ap- 
proach. While the gas algorithm uses fields, which arc discrctized on a cubic 
grid, the granulate algorithm describes particles in a continuum. A mapping 
is needed for the algorithms to interact. For the mapping a tent function F(r) 
is used: 

F(r) = f(x)f(y)f(z), f(x) = " M ^ " J' (36) 

[0 , 1 < \x/l\ , 

where / is the grid constant used for the discretization of the gas simulation. 

For the gas algorithm the porosity <pj and the granular velocity must 
be derived from the particle positions r^ and velocities v^, where i is the index 
of particle and j is the index for the grid node. The tent function distributes 
the particle properties around the particle position smoothly on the grid: 

cPj = l-^2F{Ti- Tj ), Uj = T ^-^2v i F(T i -T j ), (37) 

i 3 i 

where rj is the position of the grid point and the sum is taken over all particles. 

For the computation of the drag force on a particle the pressure drop 
VPi and the porosity 4>i at the position of the particle are needed. These can 
be obtained by a linear interpolation of the fields VP, and <f>j from the gas 
algorithm: 

4> i = Y i <t> j F(T j -T i ), VP i = 2VP J -F(r J --r i ), (38) 

3 3 

where the sum is taken over all grid points. Note that VPi is a continuous 
function of the particle position r^ . There are no discontinuities at all bound- 
aries. 

4.2 Application to Plug Conveying 

This method was applied to study plug conveying in both vertical [66] and 
horizontal [67] tubes. Plug conveying is a special case of pneumatic conveying, 
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where grains are driven through pipes by air flow. Plug conveying occurs 
when the flux of grains through the pipe is relatively high. Currently plug 
conveying is gaining importance in industry, because it causes a lower product 
degradation and pipeline erosion than dilute phase conveying. 

Unfortunately, current models [39, 65] of plug conveying disagree even on 
the prediction of such basic quantities as the pressure drop and the total mass 
flow, and these quantities have a great impact in industrial applications. One 
of the reasons for the lack of valid models is that it is difficult to study plugs 
experimentally in a detailed way. Usually experimental setups are limited to 
the measurement of the local pressure drop, the total mass flux and the veloc- 
ity of plugs. Simulational studies are handicapped by the high computational 
costs for solving the gas flow and the particle-particle interaction, and are 
therefore mostly limited to two dimensions. 

Using the above-described method, we were able to provide a detailed 
view of plugs. This approach provides access to important parameters like 
the porosity and velocity of the granulate and the shear stress on the wall at 
relatively low computational costs. Contrary to the experiments, it is possible 
to access these parameters at high spatial resolution and without influencing 
the process of transportation at all. Additional to plug profiles, characteristic 
curves of the pressure drop and the influence of simulation parameters can be 
measured. 

4.3 Results 




Fig. 10. A series of photos showing a plug moving upwards. The height of the 
shown tube is 9 cm, the frame rate is 30 Hz 



26 Harting, Hecht, Herrmann, McNamara 



In Fig. 10, we show a series of photos of plug conveying taken by Karl Som- 
mer and Gerhard Niederreiter of TU Miinchen. The particles are wax beads 
of diameter d = 1.41mm, density g s = 937kg/m 3 and a Coulomb coefficient 
of 0.21. The experimental transport channel is a vertical tube (PMMA) of 
length / = 1.01 m and of internal diameter D t — 7 mm. The air is injected at 
a constant flow rate of 2.2£/min at the bottom of the tube. As one can easily 
see, grains travel in clusters up the tube. 




Fig. 11. A series of simulation snapshots showing a plug moving upwards. The 
height of the shown tube is 12 cm, the frame rate is 100 Hz 



The simulations were carried out in a system that matched as closely as 
possible the experimental one. The same mode of transport was observed, as 
shown in Fig. 11. Not only is there a qualitative resemblance between Figs. 
10 and 11, but the simulations give the same value for the pressure drop as 
the experiments. The success of the model permitted a thorough study of the 
plugs to be carried out. For example, so-called "characteristic curves", where 
pressure drop is displayed as a function of gas velocity, could be calculated. 
The simulations also allow the study of the effects of parameters not easily 
controlled experimentally, such as the air viscosity and particle friction. The 
speed, density, size, and number of plugs were analyzed. In addition, the 
detailed structure of plugs could be studied. For example, the variation of 
density, velocity, and different components of the stress tensor were evaluated 
inside the plugs. All this information should help researchers to develop better 
models of plug conveying. 
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5 Conclusion 

In this paper we have discussed the properties of various simulation techniques 
for particles in fluids and demonstrated that there is no perfect candidate 
that is able to simulate all systems of interest and to utilize the available 
resources as efficient as possible. For each individual problem, one has to 
choose the method of choice carefully: while stochastic rotation dynamics is 
well suited to simulate systems like clay-like colloids where Brownian motion is 
important, the lattice Boltzmann method is not able to resolve the stochastic 
motion of the particles without modifications of the method. However, in 
cases where thermodynamic fluctuations are neglectably small, this approach 
is much more efficient than stochastic rotation dynamics. Like conventional 
Navier-Stokes solvers, the fluid flow can be resolved in great detail, but the 
lattice Boltzmann method is much easier to implement and to parallelize. It is 
of particular advantage if complicated boundary conditions like non-spherical 
particles or complex channel geometries come into play. The implementation 
of Navier-Stokes solvers on the other hand can be based on a long-standing and 
widespread experience with these techniques allowing to create very efficient 
solvers. In macroscopic systems like the movement of granular particles in 
air, the exact properties of the flow field are not necessary to understand 
experimentally observable parameters. Therefore, computationally much less 
demanding techniques like a coarse-grained description of the fluid should be 
applied. 
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